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ABSTRACT 


Matched-Field  Processing  (MFP)  and  Matched-Mode  Processing  (MMP) 
are  two  popialar  techniques  for  passively  localizing  an  underwater  acoustic 
emitter  in  range  and  depth.  One  major  drawback  of  these  techniques  has  been 
their  sensitivity  to  uncertainty  concerning  the  acoustic  environment.  Several 
methods  of  addressing  this  phenomenon  have  been  proposed  in  the  hterature, 
with  varying  degrees  of  success.  Achieving  high-quality  location  estimates 
remains  a  problem  except  in  simple  range-independent  experiments  or 
numerical  simulations.  In  this  study,  we  demonstrate  an  approach  for  robust, 
accurate  emitter  localization  in  a  highly  range-dependent  real  environment 
using  MMP.  The  main  factors  contributing  to  successful  localization  are:  1)  use 
of  the  high-resolution  Multiple  Signal  Classification  (MUSIC)  algorithm,  which 
performs  well  even  when  only  a  few  robust  modes  can  be  obtained  by  mode 
filtering;  and  2)  use  of  an  acoustic  propagation  model  incorporating  mode 
coupling,  which  is  able  to  generate  accurate  replica  fields  in  a  strongly  range- 
dependent  environment.  A  secondary  objective  of  the  study  was  to 
demonstrate  the  application  of  higher-order  statistical  estimation  techniques 
to  reduce  noise  effects.  Our  results  indicate  that  these  techniques  show 
unacceptable  sensitivity  to  noise-  and  model-induced  estimation  errors  and 
require  further  refinement  before  they  will  be  useful  in  the  underwater  acoustic 
localization  problem. 
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1.  INTRODUCTION 


Throughout  the  history  of  anti-submarine  warfare,  there  has  been 
great  interest  in  the  use  of  passive  acoustic  measurements  to  localize 
submarines.  Numerous  methods  of  estimating  the  Direction-Of- Arrival 
(DOA)  of  acoustic  emissions  from  a  target  of  interest  have  been  developed 
over  the  years  for  sonar  applications.  However,  these  techniques  are 
inherently  incapable  of  directly  determining  the  range  and  depth  of  a  target 
(emitter),  although  there  are  indirect  means  of  determining  range  by 
observing  target  DOA  as  a  function  of  time.  Because  knowledge  of  range  and 
depth  is  so  vital  in  military  applications,  there  has  been  considerable  interest 
in  developing  techniques  for  direct  determination  of  these  parameters.  One 
such  technique,  which  has  attracted  considerable  attention  in  recent  years,  is 
a  generalization  of  DOA  estimation  known  as  Matched-Field  Processing 
(MFP),  along  with  a  variation  on  MFP  known  as  Matched-Mode  Processing 
(MMP). 

Because  of  the  similarity  between  DOA  estimation  and  MFP/MMP, 
many  of  the  techniques  used  in  DOA  estimation  may  be  generalized  for  use  in 
MFP/MMP.  Two  of  the  most  popular  DOA  estimation  methods — ^Bartlett  and 
Minimum-Variance — ^have  been  studied  extensively  in  the  context  of  MFP 
(although  the  Minimum-Variance  method  has  not  been  addressed  in  MMP) 
[Refs.  1,  2,  3,  4,  5,  6].  The  MUSIC  method  has  received  extensive  coverage  in 
the  DOA  estimation  literature,  but  relatively  scant  attention  in  the  MFP 
literature  [Refs.  7,  8,  9]  and  no  attention  in  the  MMP  literature,  possibly 
because  of  a  perception  that  it  would  not  perform  well  in  realistic  underwater 
acoustic  environments. 
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One  of  the  most  noteworthy  results  from  MFP/MMP  studies  to  date  is 
the  great  sensitivity  of  the  location  estimates  to  uncertain  knowledge  of  the 
acoustic  environment.  Numerous  researchers  [Refs.  10,  11,  12,  13,  14,  15] 
have  studied  this  phenomenon  and  proposed  various  methods  for  addressing 
it,  with  varying  degrees  of  success.  The  problem  remains  an  open  issue. 
Because  of  this  sensitivity,  and  also  because  of  the  hmited  availability  of  real 
data  sets,  most  MFP/MMP  research  to  date  has  involved  either  simple  range- 
independent  experiments  or  numerical  simulations. 

The  objectives  of  this  dissertation  are  to: 

•  Demonstrate  the  application  of  MFP/MMP  to  experimental  data 
obtained  in  a  strongly  range-dependent  environment  during  the 
1992  Barents  Sea  Polar  Front  Experiment  [Refs.  16,  17,  18,  19]; 

•  Demonstrate  that  a  coupled-mode  propagation  model  is  able  to 
model  the  acoustic  fields  used  in  MFP/MMP  with  sufficient 
accuracy  for  high-quality  localization  estimates; 

•  Demonstrate  that  the  high  resolution  of  the  MUSIC  algorithm  in 
combination  with  MMP  is  able  to  produce  accurate  location 
estimates  in  a  realistic  environment;  and 

•  Demonstrate  the  application  of  higher-order  statistics  to  the 
MFP/MMP  problem.  Although  such  methods  have  received  some 
attention  in  the  DOA  estimation  literature  [Refs.  20,  21,  22],  they 
have  not  yet  been  applied  to  MFP/MMP. 

Chapter  II  gives  an  overview  of  pertinent  background  material, 
including  DOA  estimation  algorithms,  modeling  of  underwater  sound 
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propagation,  and  MFP/MMP  theory.  In  that  chapter,  we:  describe  three  of  the 
most  popular  DOA  estimation  algorithms — the  Bartlett,  Minimum-Variance, 
and  MUSIC  methods;  review  the  modeling  of  the  acoustic  field  via 
decomposition  into  normal  modes  for  both  range-dependent  and  range- 
independent  environments;  discuss  the  application  of  higher-order  statistics 
to  DOA  estimation;  and  derive  the  extension  of  DOA  estimation  techniques  to 
MFP/MMP.  Except  for  the  portion  regarding  application  of  the  MUSIC 
algorithm  and  higher-order  statistics  to  MMP,  this  chapter  contains  no 
original  material.  Chapter  III  gives  a  brief  overview  of  those  features  of  the 
Barents  Sea  Polar  Front  Experiment  which  are  relevant  to  this  dissertation, 
including  the  physical  characteristics  of  the  channel  and  a  description  of  the 
emitter  and  receiver.  Chapter  IV  provides  additional  information  concerning 
the  data  analysis  procedure;  it  also  presents  and  interprets  the  results  of  the 
analysis.  Chapter  V  gives  the  conclusions  reached  from  the  research  and 
proposes  areas  for  further  investigation. 
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n.  BACKGROUND 


This  chapter  provides  the  framework  for  our  anal3d;ical  approach.  It 
includes  overviews  of  the  following  topics:  DOA  estimation  fundamentals, 
including  an  extensive  description  of  the  MUSIC  algorithm;  application  of 
higher-order  statistics  to  DOA  estimation;  normal  mode  modeling  of  the 
acoustic  field;  and  theoretical  foundations  of  MFP/MMP. 

A  DOA  ESTIMATION 

One  of  the  most  fundamental  parameters  of  interest  in  military 
applications  is  the  Direction-of-Arrival  (DOA)  corresponding  to  a  target  of 
interest.  This  parameter  is  one  of  the  primary  outputs  of  virtually  all  military 
radar  and  sonar  systems.  DOA  may  be  expressed  in  terms  of  azimuth 
(bearing)  and/or  elevation.  Over  the  years,  numerous  techniques  have  been 
developed  for  DOA  estimation  (for  a  good  overview,  see  [Ref.  23]  and 
references  therein);  the  most  important  of  these  techniques  are  discussed 
briefly  in  this  chapter.  As  we  will  see  later,  these  techniques  may  be 
generalized  in  a  straightforward  manner  for  use  in  Matched-Field  Processing. 

1.  Signal  Model 

As  is  generally  the  case  in  signal  processing  algorithms,  DOA 
estimation  techniques  make  certain  assumptions  about  the  signals  being 
processed.  They  assume,  in  particular,  that  the  sound  pressure  field  in  the 
underwater  acoustic  channel  may  be  expressed  as  a  plane  wave  (i.e.,  that  the 
surfaces  of  constant  phase  are  planar).  As  we  will  see  later,  this  assumption 
is  not  true  in  general,  but  is  useful  in  many  cases  of  practical  interest.  We 
will  also  assume  that  the  signal  is  temporally  narrowband  with  center 
frequency  co;  i.e.,  that  amplitude  and  phase  modulation  do  not  introduce 
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appreciable  amplitude  and  phase  changes  over  the  physical  extent  (length)  of 
the  receiving  array.  If  a  signal  does  not  satisfy  this  latter  condition,  the 
signal  may  be  decomposed  via  Fourier  techniques  into  narrowband 
components  which  do  satisfy  the  condition.  For  the  sake  of  generality,  during 
most  of  the  background  discussions,  we  will  allow  the  number  of  emitters  to 
be  arbitrary,  even  though  the  presence  of  a  single  emitter  will  be  assumed 
dining  all  actual  data  analysis. 

For  mathematical  convenience,  we  will  conduct  our  analysis  using  the 
complex  envelope  representation  of  the  signal  rather  than  the  real  (physical) 
signal  itself.  Thus,  if  s(t)  is  the  real  signal,  the  corresponding  analytic  signal 
or  pre-envelope  [Ref.  24]  is  given  by 


where  s(^)  is  the  Hilbert  transform  of  s{t). 


The  pre-envelope  of  a  real  bandpass  signal  at  center  frequency  co  may  be 
determined  as  follows  [Ref  25]: 


•  Multiply  the  real  signal  s{t)  by  the  complex  carrier  exp(7ffl^); 

•  Pass  the  resulting  signal  through  a  low-pass  filter  to  remove  the 
component  at  twice  the  carrier  frequency; 

•  Multiply  the  resulting  baseband  signal  by  the  complex  carrier. 

All  signals  appearing  in  the  sequel  will  be  the  pre-envelopes  of  real 
narrowband  signals  unless  otherwise  indicated.  The  narrowband  assumption 
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mentioned  earlier  is  equivalent  to  requiring  the  complex  pre-envelope  of  the 
received  signal  to  be  of  the  form 

s(^)  =  Si  exp{j(ot), 

where  the  (complex)  amplitude  is  a  slowly  varying  function  of  time,  i.e., 

where  I  is  the  length  of  the  array  and  v  is  the  speed  of  soimd. 

The  concept  of  an  array  response  vector  is  one  that  arises  often  in  DOA 
estimation  and  MFP.  To  illustrate  this  concept,  we  consider  an  arbitrary 
receiving  array  of  N  elements.  We  may  represent  the  signal  as  an  iV- 
dimensional,  time-varying,  complex  vector  whose  components  are  the  signals 
at  the  individual  array  elements.  For  convenience,  we  will  assume  that  only 
the  azimuth  6  of  the  target  is  of  interest,  although  extension  to  include 
dependence  on  elevation  is  straightforward.  In  DOA  estimation,  the  array 
response  vector  a(0)  is  defined  to  be  the  unit-normalized  (i.e.,  length  of 
a(0)=l),  noise-free,  pressure  field  vector  expected  (based  on  the  modeling 
assumptions)  at  the  receive  array  given  that  an  emitter  is  at  the  angle  6.  The 
set  of  a  vectors  for  all  possible  values  of  B  is  known  as  the  array  manifold. 
More  generally,  a  could  be  a  function  of  parameters  other  than  DOA  as  well. 
For  the  simple  case  involving  a  single  emitter,  a  receiving  array  of  N 
identical,  omnidirectional  elements  in  an  arbitrary  geometry  and  a 
narrowband,  plane-wave  signal  with  center  frequency  (o,  a(0)  may  be 
expressed  as 

a(e)  =  ;^[l  e---  - 
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where  (a  function  of  6)  represents  the  time  delay  seen  by  the  incoming 
plane  wave  between  sensor  i  and  sensor  0  and  superscript  T  denotes  (non¬ 
conjugate)  matrix  transpose.  In  the  general  case  where  d  emitters  are 
present,  the  signal  model  used  in  DOA  estimation  is 

p  =  As  +  n,  (1) 

where  p  =  [pi(f)  p^^t)  •••  ^  =  ^2(^)  •"  are 

vectors  whose  elements  are  the  received  pressure  and  noise,  respectively,  at 
each  element  of  the  array;  s  =  [sl(^)  S2(^)  •••  is  the  vector  of  signals 

produced  by  the  d  emitters  (s/^)=Sjexp(/^o^),  with  a  complex  amplitude, 
because  of  the  narrowband  assumption);  A  =  [a(0i)  a(02)  •••  a(0^)]  is  a 
matrix  whose  columns  are  the  array  response  vectors  corresponding  to  the 
DOAs  of  the  d  emitters;  and  t  is  time.  Assuming  that  the  signal  and  noise  are 
uncorrelated,  the  signal-plus-noise  covariance  matrix  Rp  is  then  given  by 

Rp  =  E{pp^] 

=  AR,A«+R„,  (2) 

where  R^  =  and  R„  =  E[nn^]  are  the  signal  and  noise  covariance 

matrices,  respectively,  superscript  H  denotes  Hermitian  (conjugate) 
transpose  and  E  denotes  statistical  expectation.  An  estimate  of  Rp  derived 
from  the  measured  data  is  the  fundamental  quantity  used  in  virtually  all 
DOA  estimation  algorithms  studied  to  date. 

2.  Algorithms 

Numerous  signal  processing  algorithms  have  been  studied  in  the 
context  of  DOA  estimation.  Many  are  fairly  obvious  generalizations  of 
techniques  used  for  estimating  the  spectra  of  time  series. 
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a.  Weighted-sum  heamforming 

Weighted-sum  beamforming  (see,  for  example,  [Ref.  26])  is  the 
most  commonly  used  technique  (in  practice)  for  DOA  estimation,  and  is  used 
in  virtually  all  modem  military  radar  and  sonar  systems.  This  technique  is 
analogous  to  the  Finite  Impulse  Response  (FIR)  filters  used  in  time-series 
analysis.  The  output  bit)  of  the  beamformer  is  simply  a  linear  combination  of 
the  signals  received  by  the  elements  of  the  array, 

b{t)  =  w^p{t), 

where  w  is  the  vector  of  weights.  The  weight  vector  w  is  chosen  to  satisfy  a 
statistical  constraint  which  is  appropriate  for  the  given  situation.  The 
expected  value  B  of  the  output  power  of  the  beamformer  is  given  by 

R  =  £;[|6(f)f]  =  w"RpW.  (3) 

A  particular  value  for  w  generally  produces  high  gain  only  in  a  single  look 
direction,  i.e.,  for  a  single  value  of  6.  In  practice,  multiple  look  directions  are 
of  interest,  so  that  multiple  w  vectors  are  required  (this  approach  is 
analogous  to  the  use  of  a  “bank”  of  matched  filters  in,  for  example,  active 
sonar).  Thus,  in  general,  both  w  and  B  are  functions  of  0.  The  value  of  9  for 
which  B  is  maximized  is  taken  as  the  estimated  DOA  of  the  emitter;  this 
value  may  be  obtained  by  conventional  one-dimensional  search  techniques. 

The  Bartlett  beamformer  is  a  special  case  of  weighted-sum 
beamforming  in  which  the  weight  vectors  w  are  simply  the  array  response 
vectors  a  for  all  look  directions  of  interest.  When  the  noise  is  spatially 
homogeneous,  it  can  be  easily  shown  that  the  output  of  this  beamformer  has 
the  highest  possible  Signal-to-Noise  Ratio  (SNR)  of  any  weighted-sum 
beamformer.  This  beamformer  is  analogous  to  the  periodogram  [Ref.  271  used 
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in  time-series  analysis.  Since  for  the  Bartlett  beamformer  w=a(0),  the 
expected  value  of  the  output  power  of  the  beamformer  is 

B(e)  =  a^(0)Rpa(0)  (4) 

and  the  DOA  estimate  is  given  by 

6  =  SLTg  max  a^(0)Rpa(0). 

The  Minimum- Variance  method  (MVM),  also  (somewhat 
misleadingly)  known  as  the  “Maximum  Likelihood”  (ML)  beamformer  [Refs. 
28,  29]  selects  the  beamformer  weights  to  minimize  the  beamformer  output 
power,  subject  to  the  constraint  of  unity  gain  (i.e.,  zero  distortion)  in  the 
desired  look  direction  6.  Formally,  w  is  chosen  to 

minimize  w^RpW 
subject  to  w"a(0)  =  1. 

The  effect  of  this  minimization  is  to  produce  the  lowest  array  response  at 
directions  that  have  the  strongest  signals.  Thus,  this  method  is  useful  in 
reducing  the  effects  of  directional  noise  on  beamformer  performance.  The 
required  weights  may  be  easily  shown  (e.g.,  [Ref.  23])  (using  the  method  of 
Lagrange  multipliers)  to  be 

*  >  a«(e)R;‘a(e)- 

Substituting  into  (3)  and  simplifying  gives  a  beamformer  output  power  of 

a''(e)BXe)-  ® 

For  a  comparison  of  the  performance  of  the  Bartlett  and  MVM  techniques, 
see  Lacoss  [Ref  29]. 
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h.  MUSIC 

The  MUSIC  (for  Multiple  Signal  Classification)  method  [Refs. 
30,  31,  32]  is  one  of  the  first  and  probably  the  most  widely  studied  of  a  class 
of  techniques  which  address  the  DOA  estimation  problem  from  a  geometric 
perspective;  such  methods  offer  potentially  large  improvements  in  resolution 
with  respect  to  beamforming  [Ref.  23].  Because  MUSIC  is  central  to  the 
investigations  of  this  dissertation,  a  full  discussion  of  it  is  provided  in  the 
sequel.  This  discussion  will  attempt  to  develop  an  intuitive  understanding  of 
the  algorithm  rather  than  to  provide  a  rigorous  proof  of  its  methods. 

Signal  Subspace.  Figure  1  is  a  geometric  illustration  of  the 
behavior  of  the  received  signal  vector  p(t)  due  to  d  emitters  at  various  DOAs. 
The  coordinate  axes  represent  the  responses  of  the  N  sensors.  The 
components  of  the  vector  p  are  the  outputs  of  the  individual  sensors  and  are 
functions  of  time.  As  time  progresses,  the  tip  of  p  then  sweeps  out  a  curve  as 
shown  in  the  figure.  Due  to  obvious  limitations,  our  illustrations  can  show  at 
most  N=Z  and  will  show  sensor  outputs  as  real.  These  limitations  will  not 
adversely  affect  the  illustration  of  key  concepts. 

Sensor  3 


Sensor  2 

Sensor  1 

Figure  1:  Behavior  of  signal  vector 
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Now  consider  an  example  involving  a  single  emitter  (rf=l;  Figure 
2).  For  this  case,  as  time  progresses,  each  component  of  p  is  multiplied  by  the 
same  time-dependent  phase  factor.  The  (complex)  magnitude  of  p  changes, 
but  its  direction  does  not.  Therefore,  the  “curve”  swept  out  by  p  in  this  case  is 
a  straight  line  passing  through  the  origin. 

Sensor  3 


Sensor  2 

Sensor  1 

Figure  2:  Signal  Subspace  (c?=l) 

Next  consider  an  example  involving  two  emitters  {d=2)  at  angles 
6^  and  6^  (Figure  3).  For  this  case,  p  is  the  vector  sum  of  contributions  from 
the  individual  signals;  specifically, 

P  =  Si(0a(^i)  +  «2(0a(^2)-  (6) 

These  contributions  Sj(f)a(0i)  are  vectors  whose  magnitude 
varies  with  time,  but  whose  direction  is  fixed  by  a(0^).  Since  the  two  vectors 
are  multiplied  by  different  time-varying  phase  factors  s/0=-S,exp(/iaf) 
(provided  that  the  slowly  veirying  complex  amplitudes  are  not  perfectly 
correlated),  their  sum  p  sweeps  out  a  curve  which  is  confined  to  a  plane 
passing  through  the  origin. 
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In  general,  for  d  emitters,  the  tip  of  the  received  signal  vector  p 
sweeps  out  a  curve  which  is  confined  to  a  d-dimensional  subspace  of  <P^ .  This 
subspace  is  known  as  the  signal  subspace.  Intuitively,  then,  the  number  of 
emitters  could  be  determined  by  measuring  the  dimension  of  the  signal 
subspace.  We  will  show  later  how  this  can  be  done. 


Sensor  2 


Sensor  1 

Figure  3:  Signal  Subspace  id-2) 


Array  Manifold.  Figure  4  illustrates  an  array  manifold,  i.e., 
the  locus  of  the  array  response  vectors  a(0)  for  all  possible  values  of  6.  Note 
that,  by  definition,  a  is  independent  of  time,  so  that  time  does  not  participate 
in  this  illustration.  In  practice,  it  is  not  necessary  to  have  an  analytical 
expression  for  a(0);  we  can  instead  determine  it  experimentally  at  a  finite 
number  of  points,  store  the  results,  and  recall  them  when  desired 
(interpolating  when  necessary). 
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It  may  happen  that  the  array  manifold  “runs  over  itself’  {i.e.,  the 
mapping  a(0)  from  the  interval  0  to  2n  into  is  not  one-to-one).  In  such  a 
case,  the  array  is  said  to  have  an  ambiguity.  Such  an  ambiguity  is  not  the 
only  kind  possible.  For  example,  when  the  a  vectors  corresponding  to  3 
different  DOAs  lie  in  a  single  plane,  the  array  also  has  an  ambiguity  (for 
reasons  which  may  not  be  obvious  at  this  point).  In  general,  when  the  a 
vectors  corresponding  to  n+1  different  DOAs  lie  in  a  subspace  of  dimension  n 
or  less,  the  array  is  said  to  have  a  rank  n  ambiguity. 

Sensor  3 


Sensor  2 

Sensor  1 

Figure  4:  Array  manifold 

DOA  Determination.  For  the  case  of  one  emitter,  it  is  clear 
that  the  array  response  a  corresponding  to  the  emitter  DOA  lies  in  the  (1- 
dimensional)  signal  subspace  illustrated  in  Figure  2.  In  order  to  determine 
the  DOA,  we  would  observe  the  behavior  of  the  signal  vector  p  (if  no  noise 
were  present)  to  determine  the  signal  subspace  and  find  the  single  unit  vector 
that  spans  the  subspace.  This  unit  vector  represents  a  point  on  the  array 
manifold  corresponding  to  the  angle  of  the  emitter.  We  would  then  invert  the 
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mapping  a(0)  (which  is  one-to-one  unless  the  array  has  an  ambiguity)  to 
determine  the  DOA  6. 

Now  consider  the  case  of  two  emitters  and  recall  that  the  signal 
vector  p  is  the  sum  of  (vector)  contributions  from  the  individual  emitters  (6). 
At  any  time  t,  these  contributions  are  scalar  multiples  of  the  a  vectors 
corresponding  to  the  two  emitter  DOAs.  Thus,  the  signal  subspace  is  spanned 
by  two  vectors  &(0j)  and  a(02)  which  correspond  to  the  points  where  the  array 
manifold  intersects  the  signal  subspace  (Figure  5).  If  the  array  has  no 
ambiguities,  there  is  no  third  a  vector  that  lies  in  the  subspace.  Thus,  just  as 
in  the  case  of  one  emitter,  we  can  invert  the  a(0)  mapping  provided  by  the 
array  manifold  to  determine  the  DOAs  (Figure  5). 


Figure  5:  DOA  Determination  for  2  emitters 

In  general,  we  observe  the  data  vector  p,  determine  the  signal 
subspace,  find  its  intersection  (d  different  a  vectors)  with  the  array  manifold, 
and  invert  the  mapping  a(0)  to  estimate  the  DOAs  of  the  emitters.  It  is 
apparent  that,  in  general,  this  method  requires  the  number  of  emitters  d  to 
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be  fewer  than  the  number  of  array  elements  N  (if  d>Ar,  the  entire  array 
manifold  may  lie  within  the  signal  subspace).  Clearly,  if  the  array  manifold  is 
such  that  the  a  vectors  corresponding  to  the  emitter  DOAs  are  not  linearly 
independent  (array  ambiguity),  the  signal  subspace  dimension  is  less  than 
the  number  of  emitters.  In  such  a  case,  the  method  described  above  fails  to 
give  the  correct  number  of  emitters  and  fails  to  identify  one  or  more  of  the 
emitter  DOAs.  Such  ambiguities  can  often  be  avoided  by  proper  array  design 
(although  sometimes  physical  constraints,  such  as  in  line  arrays,  preclude 
such  design). 

Multipath.  It  is  possible  for  the  emissions  from  a  single  emitter 
to  arrive  by  two  or  more  different  paths  {e.g.,  as  a  result  of  reflection  from  the 
water  surface,  in  the  sonar  case).  To  illustrate  the  effect  of  such  a  situation, 
we  consider  the  case  of  an  array  receiving  signals  from  one  emitter  via  two 
different  DOAs.  In  such  a  case,  the  two  signals  will  exhibit  perfect  temporal 
correlation  (at  least  in  theory).  Recall  that  it  is  the  independent  variation 
with  time  of  the  array  output  due  to  individual  signals  from  different  DOAs 
that  causes  the  signal  vector  to  sweep  out  a  two-dimensional  subspace 
(Figure  3).  In  the  multipath  case,  however,  the  array  no  longer  responds 
independently  to  signals  received  from  the  two  DOAs:  the  signal  subspace  is 
one-dimensional.  In  addition,  note  that  this  signal  subspace  is  not  spanned  by 
any  combination  of  vectors  in  the  array  manifold.  Thus,  the  geometric 
approach  is  incapable  of  dealing  with  perfectly  coherent  multipath. 
Fortunately,  such  multipath  is  rarely,  if  ever,  observed  in  real  life.  However, 
the  performance  of  most  geometrically-based  algorithms  degrades  rapidly  as 
the  correlation  between  emitters  (or  multipath  arrivals)  approaches  1. 

Noise.  The  above  discussion  assumes  that  no  noise  is  present. 
With  noise  present,  the  signal-plus-noise  vector  p  sweeps  out  a  curve  that  no 
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longer  lies  in  the  signal  subspace.  Thus,  the  simple  techniques  described  are 
not  valid  if  noise  is  present. 

Essentially,  the  technique  described  above  for  noise-free 
conditions  must  be  modified  as  follows.  For  simplicity,  we  will  assume  that 
the  noise  is  spatially  isotropic  and  uncorrelated,  so  that  (2)  becomes 

R^=^AR,A^+<yX, 

where  1^^^  is  the  NxN  identity  matrix  and  the  covariance  matrices  and  Rp 
represent  theoretical,  not  estimated  statistics.  We  now  consider  the  following 
eigenvalue  problem: 

RpO  =  Ae. 

Since  Rp  is  Hermitian,  the  eigenvalues  are  real  and  the  associated 
eigenvectors  may  be  selected  to  be  orthonormal  [Ref  33].  Then, 

0  =  det[Rp-AI^] 

=  det[AR,A"  -  (A  -  <T=')I^].  (7) 

Per  the  definition  of  the  Ny.d  matrix  A,  the  columns  of  A  are  the  a  vectors 
corresponding  to  the  d  emitters.  As  mentioned  above,  if  the  array  has  no 
ambiguities,  these  a  vectors  are  linearly  independent  and  therefore  A  has  full 
(colcunn)  rank  d.  The  elements  R^Uj)  of  Rg  have  the  form  r^jS^Sj*,  where  is  a 
correlation  coefficient  between  the  ith  and  jth  signals  and  S,  represents  the 
(complex)  signal  amplitude  (not  a  function  of  time  in  this  case)  of  the  ith 
signal.  Thus,  the  dxd  matrix  R^  has  full  rank  d  unless  one  or  more  pairs  of 
emitters  are  perfectly  correlated  (in  which  case  r~l  for  some  i,  j  and 
therefore  the  ith  and  yth  columns  are  linearly  dependent).  Thus,  assuming  no 
array  ambiguities  and  less  than  perfectly  correlated  emitters,  the  first  term 
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in  equation  (7)  is  an  NxN  matrix  of  rank  d.  That  term  therefore  has  d  non¬ 
zero  eigenvalues  and  N-d  zero  eigenvalues.  Since  the  eigenvectors  of  this 
term  are  also  eigenvectors  of  the  second  term  (identity  matrix),  the 
eigenvalues  of  the  sum  of  the  two  terms  in  equation  (7)  are  the  sums  of  the 
respective  eigenvalues;  i.e.,  the  eigenvalues  of  Rp  are 


{A,-}  =  {v.+<I^ . V^+(T^<T^...,<T*}. 

The  eigendecomposition  of  Rp  may  thus  be  expressed  as 


Rp  =  EAE^ 


—  EgNgEg  +  (T  j 

where  E  is  the  matrix  of  eigenvectors,  A=diag(A,),  Eg  contains  the  columns  of 
E  corresponding  to  non-zero  eigenvalues  of  AR^A^,  and  Ng=diag(v-).  Now 
consider  the  equality 


ARsA^=  EgNgEg  (8) 

Obviously,  both  the  dxd  diagonal  matrix  Ng  and  the  Nxd  matrix  Eg  have  rank 
d.  Thus  we  see  that  the  right  hand  side  of  (8)  consists  of  a  NxN  matrix  with 
rank  d,  each  of  whose  columns  is  a  linear  combination  of  the  d  columns  of  Eg. 
Therefore,  the  N  columns  of  the  right  hand  side  must  span  the  same  subspace 
as  the  d  columns  of  Eg.  Similarly,  the  columns  of  the  left  hand  side  of  (8) 
must  span  the  same  subspace  as  the  d  columns  of  A.  Consequently,  the 
columns  of  Eg  must  span  the  signal  subspace  (the  subspace  spanned  by  the 
columns  of  A).  In  the  presence  of  noise,  then,  the  signal  subspace  may  be 
determined  by  performing  an  eigendecomposition.  The  N-d  eigenvectors 
corresponding  to  the  zero  eigenvalues  of  AR^A^  span  what  is  referred  to  as 
the  noise  subspace  (the  orthogonal  complement  to  the  signal  subspace).  In  the 
case  where  the  noise  is  not  isotropic  and  spatially  uncorrelated,  a  generalized 
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eigendecomposition  [Ref.  23]  must  be  used;  such  a  decomposition  requires 
that  an  estimate  of  the  noise  covariance  be  available.  The  existence  of  an 
estimate  of  the  noise  covariance  will  not  be  assumed  in  the  data  analysis  to 
follow,  so  we  will  not  address  this  generalized  eigendecomposition  further. 
Obviously,  there  is  some  price  to  be  paid  for  not  incorporating  the  structure  of 
the  noise  coveiriance  in  our  method;  we  expect  this  price  to  increase  as  signal- 
to-noise  ratio  decreases. 

In  the  real  world,  the  covariance  matrices  are  never  known 
exactly,  but  must  be  estimated  from  observations.  In  the  data  analysis  to  be 
presented  later,  the  following  estimate  of  the  covariance  matrix  will  be  used: 

R,  =  ??"/!,, 


where  L  is  the  number  of  observations  and  P  is  an  iVxL  matrix  whose 
columns  are  observations  of  the  received  signal  p  at  successive  times;  i.e., 

Because  of  estimation  errors,  the  N-d  noise  eigenvalues  of  this  estimated 
covariance  matrix  will  not  have  exactly  the  same  value  o^,  so  that  even  the 
number  of  emitters  d  cannot  be  estimated  with  certainty.  For  the  purpose  of 
the  present  discussion,  however,  we  assume  that  d  is  known.  Estimates  of  the 
signal  subspace  may  then  be  obtained  via  the  familiar  Linear  Least  Squares 
and  Maximum  Likelihood  techniques  [Ref.  23].  For  the  simple  case  of 
Gaussian  noise,  these  techniques  give  the  same  result,  but  one  which  is  not 
computationally  feasible  in  most  practical  situations  (since  both  involve  a 
global  minimization  over  a  space  with  dimension  equal  to  the  number  of 
emitters).  The  MUSIC  algorithm  discussed  in  the  sequel  arose  out  of  the  need 
to  reduce  computational  complexity. 
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Estimating  the  Subspace  Dimension.  In  an  optimal 
estimator,  the  problem  of  estimating  the  subspace  dimension  cannot  be 
decoupled  from  that  of  estimating  the  subspace  itself.  However,  in  the 
interest  of  reducing  computational  load,  we  can  estimate  the  dimension 
separately.  This  feature  of  the  MUSIC  algorithm  therefore  causes  it  to  be 
suboptimal. 

Estimating  the  Signal  Parameters.  As  mentioned  earlier,  in 
the  presence  of  noise,  we  can  no  longer  depend  on  precise  intersections 
between  the  signal  subspace  determined  from  the  estimated  covariance 
matrix  Rp  and  the  array  response  vectors  corresponding  to  the  signal  DOAs. 
To  estimate  these  DOAs,  we  must  therefore  determine  those  a  vectors  which 
are  “closest”  (in  some  sense)  to  the  signal  subspace.  There  are  several 
possible  methods  for  doing  this;  we  consider  only  the  simplest  method 
(Conventional  MUSIC)  here. 

Recall  from  the  earlier  discussion  that  when  the  theoretical  (i.e., 
not  estimated)  covariance  matrix  Rp  is  used  to  determine  the  signal  and  noise 
subspaces,  the  array  response  vectors  corresponding  to  the  signal  DOAs  span 
the  signal  subspace  and  are  orthogonal  to  the  noise  subspace.  Thus,  the 
squared  length  of  the  projection  of  an  a  vector  onto  the  noise  subspace,  given 
by 

Length  squared  =  a^(0)ENENa(0), 

will  be  zero  when  6  is  one  of  the  emitter  DOAs.  However,  because  of 
estimation  errors  (as  well  as  because  of  inaccuracies  in  the  signal  model  used 
to  determine  a(0)),  when  the  estimated  covariance  matrix  Rp  is  used  to 
determine  the  signal  and  noise  subspaces  (i.e.,  to  determine  E^),  this  quantity 
will  generally  not  be  zero  for  any  a.  We  must  therefore  search  the  array 
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manifold  for  the  set  of  d  array  response  vectors  which  result  in  the  lowest 
value  for  the  quantity.  An  alternative  method  for  accomplishing  this  is  to 
define  the  function 

"  a''(6/)ENE^a(6>)  ‘ 

The  estimates  of  the  signal  DOAs  then  correspond  to  the  peaks  of  this 
function.  Although  this  function  is  similar  in  some  respects  to  beamformer 
output  power  functions  such  as  that  in  (5),  the  peak  heights  of  this  function 
do  not  necessarily  provide  any  information  about  the  power  in  their 
respective  components. 

B.  HIGHER-ORDER  STATISTICS 

In  recent  years,  there  has  been  increasing  interest  in  the  use  of  higher- 
order  {i.e.,  order  greater  than  2)  statistics  in  signal  processing.  One  reason  for 
this  interest  is  that  the  statistics  known  as  cumulants  are  identically  zero  for 
Gaussian  random  variables,  provided  that  the  order  of  the  cumulant  is 
greater  than  2.  Intuitively,  we  would  expect  that,  in  situations  where  the 
noise  is  Gaussian  but  the  signal  is  not,  cumulant-based  methods  offer 
potentially  large  performance  improvements  over  conventional  methods 
based  on  second-order  statistics,  by  removing  the  noise  without  affecting  the 
signal.  A  detailed  treatment  of  higher-order  statistics  is  beyond  the  scope  of 
this  dissertation;  for  the  purposes  of  this  discussion,  a  brief  consideration  of 
the  4th-order  cumulant  will  suffice.  Broader  treatments  of  the  topic 
(including  the  material  in  the  sequel)  may  be  found  in  articles  by  Shiryaev 
[Ref.  34],  Brillinger  [Ref.  35],  and  Brillinger  and  Rosenblatt  [Ref  36];  tutorial 
articles  by  Nikias  and  Raghuveer  [Ref.  37]  and  Mendel  [Ref.  38];  and  the 
recent  textbook  by  Nikias  and  Petropulu  [Ref.  39]. 
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The  joint  characteristic  function  of  a  set  {aii,  JCg,...,  of  ti  real 
random  variables  is  defined  [Ref.  25]  by 

©2 ,  •  •  • ,  )  =  ^{exp[y(fi)i3Ci  +  (o^x^  +  •  •  •  +  «„x„)]}, 


where  E  denotes  statistical  expectation,  as  before.  The  form  of  this  function 
obtained  by  making  the  substitution  s—jfO;  is  known  as  the  moment¬ 
generating  function-,  the  moments  of  the  can  be  obtained  from  the 
coefficients  of  the  Taylor  expansion  of  the  moment-generating  function  about 
s~0.  The  second  characteristic  function  T  of  these  same  random  variables  is 
defined  as  'F=ln  O.  The  joint  cumulants  of  order  r=^j+^2+" of  these 
random  variables  are  defined  [Ref.  25]  as 


Cun^x^^  a:*- 1  -  (  ^2 . 


(10) 


a)i=a)2=..-=©n=0 


i.e.,  the  coefficients  of  the  Taylor  expansion  of 'P  about  <»~0.  For  the  4th-order 
cumulant  of  zero-mean  real  random  variables  x^ ,  it  can  be  shown  [Refs.  34, 
39]  that  (10)  reduces  to 

Cum(xi,x2,x^,xf)  =  Elx-^x^x^x^]  -  jB[xiX2]  •  E[oc^x^ 

-.^[ACiAJg]  •  ^[31:20:4]  -  Elx-^x^]  •  Elx^x^] ,  (11) 


It  may  further  be  shown  [Ref.  25]  that,  if  the  Xi  are  jointly  Gaussian,  the  4th 
order  moment  is  given  by 

•®[^1^2^3^4]  ~  -^[^1^2]  '  ■^[^3^4]  •^[^1^3]  ‘  ■^[^2^4]'*’  -^[^1^4]  '  •®[^2^3]>  (12) 

therefore  the  cumulant  vanishes  as  claimed.  For  complex  random  variables, 
(11)  takes  the  form  [Ref.  40] 
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(13) 


Cum{xy, xl,  JCg,  X4)  =  Elx^x^x^x^]  -  E[x^xl]  •  E[x^xl] 

-jB[a:ia:3]  •  E[x\x\]  -  ^^1X4]  •  Elx^x^]  ^ 

where  the  third  term  is  generally  assumed  to  be  zero  due  to  the  S5rmmetry 
property  between  the  real  and  imaginary  parts  of  a  stationary,  bandpass, 
complex  process  [Ref.  41].  We  will  retain  this  term  for  generality. 


A  spatial  4th-order  cumulant  matrix  C4  may  be  defined  as  follows: 


P2{t)pl{t) 

(14) 


where  denotes  complex  conjugate  and  the  p-{t)  are  sensor  signals  (complex 
pre-envelopes)  at  each  of  N  sensors  in  an  array.  This  definition  differs  from 
the  one  used  by  Nikias  and  Petropolu  [Ref.  39],  but  was  selected  so  that  C4  is 
Hermitian.  By  substituting  (11)  into  (14),  and  using  vector  notation,  we 
obtain 

C4  =  ■e{(p  °  p*  )(p  °  pT  }  “  -®{p  °  P’  }^{(P  °  P*  )^  } 

-.e{pp"}o£;{pV}-e{pp^}oJ5;{p*p»}  ^  ^5^ 


where  o  denotes  the  element-by-element  product  of  the  vectors  and 
superscript  T  indicates  (non-conjugate)  matrix  transpose.  The  fourth  term  of 
(15)  is  generally  assumed  to  be  zero  due  to  the  symmetry  property  discussed 
in  conjunction  with  (13)  but  will  be  retained  for  generality.  Substituting  the 
signal  model  of  equation  (1)  (for  the  case  of  a  single  emitter)  into  (15)  gives 

C4  =  o  a*s*)(as  o  a*s*)^|  -  £?|as  o  aV|.B|(as  o  a*s*)^| 

-£;{ass*a"}  o  £;{a*s*sa^}  -  E{assa^}  o  £;{a*s*s*a"} 
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provided  that  the  signal  and  noise  are  independent  (the  single-emitter 
assumption  is  made  here  for  notational  convenience  and  will  also  be  made 
during  the  actual  data  analysis  in  the  sequel).  Since  a  is  deterministic,  it  may 
be  factored  out  of  the  expectation  operators,  and  we  obtain,  after  some 
algebra 

C4=r(aoa*)(aoa*f^  (16) 

where 

7  =  Cum{s(^),  s*(^),  s{t),  s*(f)} 

is  the  kurtosis  measure  of  the  (single)  emitter  signal  s(t)  and  a  is  the  array 
response  vector  corresponding  to  the  location  of  the  emitter. 

The  structure  of  C4  is  therefore  identical  to  that  of  the  covariance 
matrix  defined  in  (2),  except  that:  1)  the  noise  covariance  vanishes  due  to  the 
properties  of  the  4th-order  cumulant;  2)  the  array  response  vector  a  is 
replaced  by  aoa*,  which  is  real;  and  3)  C4  has  rank  1  (see  (16))  due  to  the 
assumption  of  a  single  emitter  and  is  real.  This  structure  allows  the  use  of  a 
modified  version  of  the  previously  discussed  MUSIC  algorithm  as  follows: 

•  Form  an  estimate  C4  of  the  C4  matrix  from  the  measured  data 
using  (15),  with  expectation  operators  replaced  by  time  averages; 

•  Perform  an  eigendecomposition  of  C4.  Because  it  has  rank  1,  there 
will  be  (theoretically)  a  single  non-zero  eigenvalue. 

•  Select  the  eigenvector  corresponding  to  the  largest  eigenvalue;  as 
shown  eeirlier  (8),  this  eigenvector  must  span  the  same  subspace  as 
the  a  vector  corresponding  to  the  actual  emitter  location.  The 
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remaining  eigenvectors  span  the  noise  subspace  and  form  the 
columns  of  E^. 

•  Form  the  function 

"  [a(e)»a*(<E„ES[a(9).a»(9)]' 

the  estimate  for  the  emitter  location  corresponds  to  the  peak  of  this  function. 
This  method  can  be  easily  generalized  to  the  multiple-emitter  scenario, 
provided  that  the  signal  subspace  dimension  is  selected  to  correspond  to  the 
number  of  emitters.  However,  oiu"  experimental  work  with  real  data  does  not 
require  this  generalization. 


Cumulant-based  versions  of  the  Bartlett  and  MVM  methods  also  exist 
[Ref.  39],  but  these  will  not  be  used  in  the  data  analysis.  It  should  also  be 
noted  that  the  form  of  the  4th-order  cmnulant  appearing  in  (15)  is  a  reduced 
form  of  that  used  by  Porat  and  Friedlander  [Ref.  20]. 


C.  ACOUSTICS  AND  MODELING 
1.  Helmholtz  Equation 

Whereas  in  DOA  estimation  the  received  signals  are  assumed  to  be 
plane  waves,  in  MFP  the  pressure  field  amplitude  p{r,z)  in  an  underwater 
channel  (see  Figure  6)  due  to  a  narrowband  emitter  with  center  frequency  (o 
is  assumed  to  satisfy  the  Helmholtz  equation  [Ref.  42] 


1^ 

r  dr 


dp{r,z) 

dr 


(18) 


where  r  is  the  range  from  the  emitter  to  the  observation  point,  z  is  the  depth, 
and  k=(olc  is  the  wavenumber.  Cylindrical  coordinates  are  used  in  (18) 
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because  of  the  S3nnmetry  which  exists  when  changes  in  sound  speed  in  the 
azimuthal  direction  are  negligible.  The  sequel  will  present  an  overview  of 
methods  for  solving  (18)  numerically  which  are  pertinent  to  this  dissertation. 
In  this  section,  it  will  be  assumed  that  no  observation  noise  is  present. 


ro 
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Figure  6:  Generic  Underwater  Channel 


2.  Normal  Mode  Solution 

In  a  channel  where  all  properties  (sound  speed,  water  depth,  bottom 
type,  etc.)  are  independent  of  horizontal  range,  it  is  well  known  [Ref.  42]  that 
the  pressure  field  p  at  range  r  (relative  to  an  arbitrary  origin;  see  Figure  6) 
sufficiently  far  away  from  the  emitter  and  at  depth  z  may  be  expanded  in 
terms  of  normal  modes  as  follows: 


p{r,z,t;ro,Zo)='^ 


m=l 


(19) 


where:  Zq  and  are  the  emitter  depth  and  range,  respectively;  r  is  taken  to 
be  greater  than  r^;  A  is  a  constant  which  depends  on  the  emitter  power;  and 
the  and  are  the  eigenfunctions  and  eigenvalues,  respectively,  of  the 
ordinary  differential  equation 
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The  boundary  conditions  for  (20)  depend  on  the  acoustic  properties  of  the 
surface  (which  is  always  assumed  to  be  pressure-release;  i.e.,  p(z=0)=0)  and 
the  bottom.  Attenuation  due  to  sediment  and  water  column  is  incorporated, 
as  is  customary,  via  a  small  imaginary  part  in  the  k^.  Although  the  sum  in 
(19)  is  over  values  of  m  from  1  to  infinity,  all  modes  for  which  m  is  greater 
than  some  integer  M  are  attenuated  enough  to  be  ignored  at  the  ranges  of 
interest;  such  modes  include  the  strongly  bottom-interacting  and  evanescent 
modes.  In  the  discussions  to  follow,  we  will  ignore  these  modes  and 
incorporate  only  the  lowest  M  modes  in  our  normal  mode  expressions.  It 
should  be  noted  that,  in  general  (i.e.,  for  range  near  zero),  the  expansion  (19) 
is  only  approximate,  since  it  only  accounts  for  the  discrete  spectrum  of  the 
modal  solution  to  the  Helmholtz  equation  (18).  This  fact  does  not  present  a 
problem,  however,  because  the  contribution  to  the  pressure  field  from  the 
continuous  spectrum  is  negligible  at  the  ranges  of  interest  in  MFP/MMP. 
Numerous  methods  exist  for  implementing  (19)  on  a  computer  (see,  for 
example,  [Ref.  43]).  As  is  characteristic  of  solutions  of  boundary  value 
problems  such  as  (20),  the  are  orthogonal  with  respect  to  the  density 
function p(z),  i.e., 

Jo"  ^  =  vj{m  -  n),  (21) 


where 


This  orthogonality  will  prove  useful  later  in  the  discussion  of  Matched-Mode 
Processing. 

3.  Adiabatic  Approximation 

As  mentioned  above,  the  normal  mode  expression  (19)  is  strictly  valid 
only  in  range-independent  environments;  nevertheless,  it  can  be  modified 
slightly  to  apply  to  a  limited  class  of  range-dependent  environments. 
Whenever  range  dependence  exists  (due  to  a  sloping  bottom  or  change  in 
sound  speed  profile,  for  example),  it  is  clear  that  the  and  must  be 
functions  of  range.  If  the  range  dependence  is  relatively  weak,  it  can  be 
assumed  that  a  mode  does  not  exchange  energy  with  other  modes,  but  merely 
adapts  itself  to  local  environmental  conditions.  This  assumption  is  known  as 
the  “adiabatic”  approximation.  In  such  a  case,  (20)  is  solved  at  each  of  the 
ranges  of  interest.  The  resulting  and  are  then  functions  of  range,  so 
that  the  pressure  field  can  be  expressed  as 

p(r,  2,  To,  Zo)  =  S  rr'-( . ^ . n  ^o)  exp 

4  Coupled  Mode  Model 

Due  to  the  energy  exchange  between  modes  in  a  strongly  range- 
dependent  environment,  the  use  of  a  normal-mode  model  in  such  an 
environment  becomes  more  complex.  For  simplicity,  we  assume  that 
horizontal  refraction  and  azimuthal  scattering  are  negligible  (see  [Ref.  44]  for 
a  fully  three-dimensional  treatment).  As  will  be  discussed  in  more  detail 
later,  this  approximation  is  satisfactory  for  the  acoustic  environment 
addressed  in  the  data  analysis  to  follow  [Ref.  19]. 

Mode  coupling  may  be  accounted  for  in  a  manner  consistent  with  the 
notation  used  above  by  defining  a  range-dependent  mode  amplitude  K(r) 


r 

j(Ot-jjk^{r)dr 


ro 


.(22) 
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(j.e.,  the  complex  constant  A  in  (22)  becomes  a  function  of  both  range  and 
mode  number  m)  [Ref.  45].  The  pressure  field  is  then  given  by 


p(r, z, t; To, Zo)  =  S ^o) exp 


jG)t-jjk^{r)dr 

^0 


.(23) 


In  this  expression,  all  coupling  between  modes  is  accounted  for  by  the  (i.e., 
each  A^  depends  on  the  mode  amplitudes  and  phases  of  all  other  modes  along 
the  propagation  path).  The  Broadband  Coupled  Mode  model  developed  by 
Chiu  et  al.  [Ref.  19]  has  demonstrated  high  accuracy  in  predicting  the  modal 
structure  observed  in  the  Barents  Sea  Polar  Front  Experiment  and  will  be 
used  for  all  data  analysis  in  the  sequel. 

D.  MATCHED-FIELD  AND  MATCHED-MODE  PROCESSING 


A  good  overview  of  this  topic  may  be  found  in  the  textbook  by  Tolstoy 
[Ref.  46].  Only  background  material  pertinent  to  the  later  data  analysis 
(along  with  some  preliminary  material)  will  be  presented  here. 

1.  Motivation 

Several  limitations  in  appl3dng  DOA  estimation  techniques  to  the 
underwater  acoustic  localization  problem  arise  from  the  inability  of  the 
simple  plane-wave  signal  model  to  describe  the  acoustic  field  adequately.  In 
this  section,  we  describe  these  limitations,  which  provide  the  motivation  for 
study  of  MFP/MMP. 

a.  Ability  to  determine  target  parameters 
Because  of  the  assumption  that  the  received  signal  is  a  plane 
wave,  the  only  degree  of  freedom  in  the  array  response  vector  is  DOA. 
Therefore,  DOA  estimation  is  inherently  incapable  of  providing  any  other 
information.  In  particular,  it  gives  no  range  information,  which  in  military 
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applications  is  of  vital  interest  (although  it  should  be  pointed  out  that  it  is 
possible  to  determine  range  by  observing  DOA  information  over  time, 
provided  that  certain  assumptions  concerning  target  motion  are  veilid;  see 
[Ref  471). 

h.  Estimation  error 

Since  the  speed  of  sound  in  seawater  is  a  function  of  position, 
the  sound  “rays”  {i.e.,  paths  normal  to  the  surfaces  of  constant  phase)  are 
curved.  Consequently,  the  DOA  estimates  will,  in  general,  be  different  from 
the  actual  directions  of  the  emitters.  In  many  practical  situations,  the 
receiving  array  has  no  depth  extent  (e.g.,  a  horizontal  line  array)  and  the 
ocean  may  be  considered  to  be  horizontally  stratified  (i.e.,  speed  of 
propagation  is  a  function  primarily  of  depth).  In  such  situations,  the  plane 
wave  signal  model  is  relatively  accurate.  Even  in  such  cases,  the  DOA 
estimates  will  often  be  significantly  in  error  due  to  reflection  of  sound  energy 
from  boundaries  (surface  and  bottom). 

c.  Loss  of  gain 

We  have  noted  earlier  that  the  Bartlett  beamformer  results  in 
the  highest  output  SNR,  provided  that  the  noise  is  spatially  homogeneous. 
Essentially,  this  feature  is  due  to  the  fact  that  the  signals  from  different 
receive  array  elements  add  constructively  (because  the  beamformer 
introduces  phase  shifts  to  account  for  the  shape  of  the  wavefront)  while  the 
noise  components  do  not.  In  an  underwater  acoustic  channel,  however,  the 
wavefronts  may  not  be  planar  (particularly  in  the  vertical  direction),  so  that 
the  signals  from  the  array  elements  will  no  longer  be  added  perfectly 
constructively  (since  the  beamformer  assumes  the  wavefront  is  planar  when 
it  actually  is  not).  Consequently,  the  gain  of  the  beamformer  will  be  lower 
than  when  the  incoming  signal  is  a  plane  wave. 
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d.  Inability  to  resolve  multipath  arrivals 

As  mentioned  earlier,  many  eigenstructure-based  methods 
(including  MUSIC)  are  unable  to  handle  signals  from  multiple  highly 
correlated  emitters.  Reflections  from  boundaries  in  the  underwater  acoustic 
problem  can  often  result  in  the  same  signal  arriving  at  the  receiver  from 
multiple  directions.  This  situation  is  equivalent  to  the  presence  of  multiple, 
highly  correlated  emitters;  thus,  MUSIC  and  similar  methods  break  down  in 
such  a  situation. 

2.  Generalization  of  array  processing  algorithms 
a.  Matched-Field  Processing  (MFP) 

For  the  case  of  a  single  signal,  (1)  becomes 

p  =  as(^)  +  n.  (24) 


We  can  use  (23)  to  express  the  pressure  at  a  receiving  array  consisting  of  a 
set  of  iV  vertically-aligned  hydrophones  at  depths  [z^,  Z2,...,Zjf}  in  the  form  (24) 
if  we  identify 


p{r,z^,t) 

p{r,Z2,t) 

p{r,Zj^,t) 


(25) 


and 


a 


-  Kir) 


Z„»(r)Z,„(ro,2o)exp 


(26) 


where  exp(7faf)=s(0  and  ZJr)=[ZJr,z^),  ZJr,Z2),...,  Z„(r,2^)]^.  The  array 
response  vector  a  is  now  a  function  of  emitter  range  and  depth  2^  rather 
than  azimuth,  as  in  DOA  estimation.  The  a(ro,2o)  are  obtained  by  calculating 
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the  pressure  field  seen  at  the  receive  array  using  an  appropriate  propagation 
model,  for  every  possible  {r^z^)  combination  of  interest.  In  the  previous 
discussions  on  DOA  estimation,  no  assumptions  were  made  concerning  the 
form  of  a.  Furthermore,  the  dependence  of  a  on  azimuth  alone  was  for  the 
purposes  of  illustration  only,  and  is  not  required  for  the  DOA  estimation 
techniques  to  be  valid.  Therefore,  the  algorithms  discussed  above  in  the 
context  of  DOA  estimation  may  be  used  in  MFP  as  well  [Refs.  3,  7].  It  should 
be  noted  that  MFP  does  not  require  that  the  acoustic  field  be  expressed  in 
terms  of  normal  modes;  the  above  analysis  is  valid  for  other  types  of 
propagation  models  as  well. 

The  well-known  principle  of  acoustic  reciprocity  (see,  e.g. ,  [Ref. 
48])  is  very  useful  in  minimizing  the  computations  required  to  generate  the 
array  manifold.  This  principle  states  that,  under  a  set  of  reasonable 
assumptions,  the  acoustic  pressure  at  a  location  (r^z)  generated  by  a  simple 
source  at  location  (r^yZg)  is  the  same  as  the  acoustic  pressure  at  location 
generated  by  that  same  source  at  location  {r,z).  In  MFP,  the  construction  of 
the  array  manifold  requires  that  the  field  at  a  known  receiver  location  (r,z)  be 
computed  for  every  possible  emitter  location  whereas  one  run  of  a 

propagation  model  will  generally  produce  the  acoustic  field  at  every  possible 
receiver  location  due  to  an  emitter  at  a  known,  fixed  location.  Thus,  it  may  be 
seen  that  this  principle  of  acoustic  reciprocity  allows  construction  of  one 
component  (corresponding  to  one  element  of  the  receive  array)  of  the  array 
manifold  with  a  single  run  of  a  propagation  model.  The  total  niunber  of  runs 
required  will  be  the  same  as  the  number  of  elements  in  the  receive  array. 
Clearly,  this  approach  allows  huge  savings  in  computation  compared  to  a 
“brute  force”  approach  which  performs  one  run  of  a  propagation  model  for 
each  possible  source  location. 
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The  procedure  for  MFP  may  be  outlined  as  follows: 


•  Determine  the  pre-envelopes  of  the  received  (real)  signals  at  each 
receive  hydrophone.  Use  these  pre-envelopes  to  generate  an 
estimate  (using  time  averages  instead  of  statistical  expectations)  of 
the  covariance  matrix  Rp  (or  the  4th  order  cumulant  matrix  C4); 

•  Using  a  suitable  propagation  model  and  invoking  the  principle  of 
acoustic  reciprocity,  generate  the  array  manifold  vectors 
(which  in  the  context  of  MFP/MMP  are  generally  called  replica 
fields)  for  values  of  (r^^)  on  a  suitable  grid; 


•  Generate  the  functions  in  equations  (4)  (Bartlett),  (5)  (MVM),  (9) 
(MUSIC),  or  ( 17)  (cvunulant  MUSIC),  as  desired.  These  will  now  be 
functions  of  two  variables  (r^yZo)  rather  than  one  (0),  i.e. , 


■8(^0,  Zq)  =  Zo)  (Bartlett), 

1 


%o>-2^o)  = 


a"(ro,Zo)R;^ro,Zo) 
1 


(27) 

(28) 


PM[ro,Zo)  = 


a  (^0 »  Zq  )Ej4E^a(ro ,  ) 

1 


(MVM), 

(MUSIC),  (29) 

(Ciunulant).  (30) 


[a(^o,  ^o)  °  a*(^o>  ^o)rENE"[a(ro,  2:0)  °  ^(r,,  z,)] 

The  estimated  emitter  locations  are  the  combinations  for 

which  these  functions  are  maximized.  In  MFP/MMP,  plots  of  these 
functions  are  generally  referred  to  as  ambiguity  surfaces, 
h.  Matched-Mode  Processing  (MMP) 

MMP  [Refs.  4,  49,  50,  51]  requires  the  pressure  field  in  the 
channel  to  be  expanded  in  terms  of  normal  modes.  As  will  become  clearer 
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later,  it  may  be  understood  intuitively  as  a  transformation  of  the 
observations  from  the  space  of  individual  hydrophones  to  the  space  of  modal 
amplitudes.  This  transformation  is  known  as  mode  filtering.  As  before,  we 
will  assume  that  an  accurate  representation  of  the  field  (23)  requires  only  a 
finite  munber  of  terms  M. 

As  with  MFP,  the  pressure  field  is  sampled  with  an  AT-element 
vertical  array  of  identical  elements  at  depths  {^l,  located  a  distance  r 

relative  to  an  arbitrary  origin.  The  vector  of  received  pressures  (defined  in 
(23)  and  (25)),  without  additive  noise,  may  be  expressed  in  matrix  form  as 


p=Zu, 


.(31) 


where 


Z  = 


Z2(r,z,) 

••• 

Z,(r,Z2) 

^2(^22) 

1 - 

... 

u  =  [mi 

U2  ••• 

,  with 

An(r) 

It  /  \ / 

-rZJro, 

Zo)exp 

jciJt-j 

and 


'b 


(32) 


Z  thus  contains  all  information  about  the  receiving  array,  while  u  contains 
all  information  about  the  location  of  the  emitter  relative  to  the  receiver.  Both 
Z  and  u  contain  information  about  the  channel  via  the  eigenfunctions 
which  are  derived  from  normal  mode  analysis.  An  estimate  u  of  u  may  be 
obtained  from  (31)  using  either  of  two  classes  of  methods. 

The  first  class  of  methods  regards  estimation  of  u  from  a  purely 
mathematical  perspective,  namely,  as  a  least-squares  problem  [Ref.  33], 
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either  overdetermined  (Ar>M)  or  underdetermined  {N<M).  In  this  problem,  u 
is  selected  to  minimize  the  quantity 


Zu- 


(33) 


the  squared  Euclidean  norm  of  the  residual  (Zu-p).  In  the  overdetermined 
case,  assuming  that  Z  has  full  rank,  the  solution  u  is  unique.  In  that  case, 
the  selection  of  a  suitable  method  is  based  on  considerations  of  numerical 
stability  and  computational  complexity.  If  the  problem  is  underdetermined 
(as  is  the  case  with  the  data  to  be  analyzed  in  this  dissertation)  or  Z  is  rank- 
deficient,  there  is  an  infinitude  of  vectors  u  which  minimize  (33).  Two 
subclasses  of  least-squares  methods  exist  in  this  case:  those  which  produce  a 
solution  u  with  minimum  norm  (such  as  the  pseudo-inverse  method  discussed 
below)  and  those  which  do  not  (such  as  certain  versions  of  the  QR 
factorization  method).  Methods  in  the  latter  subclass  give  solutions  with 
significantly  greater  sensitivity  when  p  is  contaminated  by  observation  noise. 
We  will  consider  only  the  pseudo-inverse  method  in  this  dissertation. 

The  singular  value  decomposition  (SVD)  of  Z  is  given  by  [Ref. 
33] 


U"ZV  =  diag((Ti,  <T2,  . . . ,  (Tp ) 


1 

1 

1 

1  0 

1 _ 

where  U  and  V  are  N'xN  and  MxM  unitary  matrices  respectively, 
p=min(M,iV),  and  the  are  the  (real  and  non-negative)  singular  values  (some 
of  the  CTj  will  be  zero  if  Z  is  rank-deficient).  Note  that  the  matrix  partition 
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shown  assumes  i\r<M ;  the  partition  when  N>M  is  analogous.  The 
pseudoinverse  Z*  of  Z  is  defined  as 


where 


=  diag 


(34) 


and  i2=rank(Z).  The  minimimi-norm  solution  of  (31)  can  be  expressed  as 


u=Z^p.  (35) 

Obviously,  the  expression  (30)  is  very  sensitive  to  the  presence  of  small  but 
non-zero  singular  values  (due,  e.g.,  to  roundoff  error).  This  phenomenon  may 
be  satisfactorily  dealt  with  by  treating  all  singular  values  on  the  order  of 
machine  precision  (or  smaller)  as  zero.  As  is  well  known,  if  N>M  and  Z  has 
full  rank,  the  unique  least-squares  solution  is  given  by 

u  =  Z^p  =  (Z'^Z)'Vp. 

This  method  of  modal  decomposition  is  referred  to  by  Yang  [Ref.  52]  as  the 
“Eigenvector  Method”.  Obviously,  modes  which  are  so  poorly  sampled  by  the 
receive  array  that  the  corresponding  columns  of  Z  are  nearly  linearly 
dependent  cannot  be  resolved  using  this  approach. 

As  mentioned  earlier,  the  analysis  to  follow  does  not  assume  the 
existence  of  an  estimate  of  the  noise  covariance.  It  is  worth  noting  that,  when 
p  is  contaminated  by  observation  noise  with  known  covariance  (i.e.,  known, 
except  possibly  for  a  constant  multiplicative  factor),  u  should  be  determined 
via  the  generalized  least-squares  method  [Ref.  33]  in  order  to  ensure  that 
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undue  weight  is  not  given  to  data  from  hydrophones  with  high  levels  of 
observation  noise.  Since  this  method  will  not  be  used  in  our  analysis,  it  will 
not  be  discussed  further. 

The  second  class  of  methods  (see,  for  example,  [Ref.  51])  for 
estimating  u  regards  the  problem  from  a  more  physical  perspective;  i.e.,  it 
exploits  the  orthogonality  property  of  the  mode  functions  (columns  of  Z)  (21). 
In  (21),  if  the  spacing  of  hydrophones  is  sufficiently  dense,  the  integral  may 
be  replaced  by  a  sum  without  significant  loss  of  accuracy.  Furthermore, 
within  the  water  column,  the  density  p  is  approximately  independent  of 
depth.  Thus,  the  columns  of  Z  eue  approximately  orthogonal,  at  least  for  the 
lower  modes  (i.e.,  for  those  modes  which  are  well  sampled  by  the  receiver 
depths  2;),  i.e., 

Z''Z.^diag(n,v, . v„), 

where  Az  is  the  spacing  of  the  vertical  grid  on  which  the  mode  function  is 
evaluated.  To  obtain  an  estimate  of  u,  we  premultiply  (31)  by  Z® 

Z"p  =  Z"Zu 

“;^diag(vi,V2,...V3,)-u. 


We  may  thus  take  the  estimate  of  u  to  be 


-  Az  j.  I 
u  =  — diag 

P 


V^l’^2 


Z"p. 


''M  J 


(36) 


This  method  of  estimating  u,  which  we  will  refer  to  as  the  projection  method, 
is  obviously  very  similar  mathematically  to  the  pseudoinverse  method. 
Because  of  the  mode  orthogonality  assumption,  only  modal  amplitudes 
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corresponding  to  those  columns  of  Z  which  form  a  (nearly)  mutually 
orthogonal  set  may  be  accurately  estimated. 

When  observation  noise  is  present,  (31)  becomes 

p  =  Zu+n. 

As  mentioned  earlier,  because  we  do  not  assume  that  an  estimate  for  the 
noise  covariance  is  available,  the  presence  of  this  observation  noise  does  not 
affect  how  we  estimate  u  (although,  of  course,  the  value  of  the  estimate  will 
be  affected).  Using  the  pseudoinverse  method  to  estimate  u,  we  obtain 

u  =  Z"p 
=  Z^Zu  +  Z^n 

=  u  +  Z^n  .  (37) 


From  (32)  we  have 


u  =  as(^) , 


where 


= 

iti 


a(ro,2o)=[«i(^o^o). 
_  A»(^;^o>^o) 


^K{r){r-r^ 


^™(^o.2o)exp 


-j\K(r)dT 


,  and 


(38) 


s(^)=exp(7<a^). 


If  we  define 


n'  =  Z  n. 
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(37)  becomes 


u  =  as(^)  +  n'.  (39) 

The  expression  (39)  is  of  the  form  (1),  again  for  the  special  case  of  a  single 
signal.  Therefore,  once  u  is  known,  the  MFP  techniques  discussed  above  may 
be  used  to  find  the  emitter  range  and  depth.  As  discussed  earlier,  no  estimate 
of  the  observation  noise  will  be  used  in  the  analysis  to  follow.  Thus,  the  fact 
that  the  noise”  n  has  a  different  character  from  n  does  not  affect  our 
analysis.  The  analysis  for  the  case  of  the  projection  method  is  essentially  the 
same  as  the  foregoing  and  will  not  be  presented  separately.  Yang  [Ref.  52] 
notes  that  this  method  of  modal  decomposition  gives  a  localization  estimate 
which  is  mathematically  equivalent  to  the  MFP  approach  when  the  Bartlett 
processor  and  all  modes  are  used. 

Regardless  of  which  method  of  estimating  u  is  used,  care  is 
required  in  selecting  which  subset  (i.e.,  which  components  of  u  and  a)  of  the 
full  mode  set  (obtained  from  either  (35)  or  (36))  is  to  be  used,  because: 

•  As  the  mode  number  increases,  so  does  the  vertical  wavenumber 
[Ref.  42];  thus,  for  higher-order  modes,  a  closer  vertical  spacing  of 
receive  hydrophones  is  necessary  for  adequate  “sampling” 
(analogous  to  the  sampling  theorem  in  time-series  analysis); 

•  Propagation  models  generally  show  greater  sensitivity  to 
uncertainties  in  the  environmental  parameters  when  predicting 
high-order  modes  than  when  predicting  low-order  modes  (see,  e.g., 
[Ref  521); 
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•  Modes  which  are  only  weakly  present  {i.e.,  which  have  low  modal 
amplitude)  at  the  receiver  location  can  generally  not  be  estimated 
accurately; 

•  Only  those  modes  which  have  most  of  their  energy  at  depths  within 
the  physical  extent  of  the  receiving  array  are  likely  to  be  estimated 
accurately  (see,  e.g.,  [Ref.  52]); 

•  When  the  number  of  receive  hydrophones  is  less  than  the  number  of 
modes  supported  by  the  channel  (as  is  the  case  with  the  data  to  be 
analyzed  later  in  this  dissertation),  the  inversion  of  (31)  is  an 
underdetermined  problem  and  therefore  cannot  provide  accurate 
values  for  all  modes. 

The  first  and  second  considerations  usually  favor  the  lower-order  modes.  This 
generalization  is  not  valid  in  all  situations:  for  example,  incorporation  of 
poorly  resolved  higher-order  modes  into  the  estimator  can  sometimes  reduce 
sidelobe  levels  when  the  low-resolution  Bartlett  processor  is  used  [Ref.  52]. 
The  third  consideration  also  tends  to  favor  low-order  modes,  at  least  when 
the  emitter  and  receiver  are  widely  separated  (since  higher-order  modes  are 
attenuated  more  rapidly).  The  fourth  consideration  is  relatively  easy  to 
employ,  since  the  modal  structure  (i.e.,  mode  shapes)  at  the  receiver  location 
is  known.  The  fifth  consideration  favors  modes  which  are  well  sampled  by  the 
receiving  array  and  which  are  therefore  nearly  orthogonal,  since  only  nearly 
orthogonal  modes  are  well  resolved  by  mode  filtering.  Again,  these  modes  are 
generally  the  low-order  ones.  In  general,  an  initial  localization  estimator 
should  be  constructed  based  on  a  mode  set  selected  using  the  above 
considerations.  The  peaks  of  this  estimator  may  be  regarded  as  candidate 
emitter  locations.  Then,  revised  estimators  (one  estimator  per  candidate 
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source  location)  may  be  obtained  by  using  only  those  modes  which  are 
expected  to  be  present  at  the  receiver  due  to  sources  at  these  candidate 
locations.  Thus,  a  strategy  of  iterative  improvement  may  be  used  to  refine  the 
estimator. 

The  MMP  technique  may  be  summarized  as  follows: 

•  Perform  quadrature  demodulation  on  the  received  (real)  signals  at 
each  receive  hydrophone  to  obtain  p(f); 

•  Use  either  the  pseudo-inverse  or  projection  methods  to  estimate  u(^) 
fi*om  p(f)  using  (31); 

•  Generate  an  estimate  of  the  modal  covariance  matrix  R^=S[uu®] 
(or  the  4th  order  modal  ciimulant  matrix  C4); 

•  Using  a  suitable  propagation  model  and  invoking  the  principle  of 
acoustic  reciprocity,  generate  the  array  manifold  vectors 
(from  (38))  for  values  of  on  a  suitable  grid; 

•  Select  a  subset  of  the  full  mode  set  for  use  in  further  processing; 

•  Generate  the  functions  in  equations  (27)  (Bartlett),  (28)  (MVM), 

(29)  (MUSIC),  or  (30)  (cumulant  MUSIC),  as  desired.  The  estimated 
emitter  location  is  the  combination  for  which  the  function  is 

maximized. 

•  As  is  apparent  from  our  discussion  above  concerning  mode 
selection,  the  major  advantage  (at  least  from  the  perspective  of  our 
research)  of  performing  MFP  in  mode  space  (i.e.,  MMP)  rather  than 
in  hydrophone  space  is  that  estimation  errors  due  to  environmental 
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mismatch  may  be  reduced  by  using  by  using  only  robust  modes, 
which  are  less  sensitive  to  such  mismatch. 
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in.  EXPERIMENTAL  SETUP 


This  chapter  provides  an  overview  of  those  aspects  of  the  Barents  Sea 
Polar  Front  Experiment  which  are  applicable  to  this  research.  This 
experiment  provided  the  data  on  which  the  MFP/MMP  algorithms  were 
tested. 

A  ENVIRONMENT 

The  data  used  in  the  analysis  to  follow  was  obtained  during  the  1992 
Barents  Sea  Polar  Front  Experiment  [Refs.  16, 17, 18, 19].  Most  of  the  details 
of  that  experiment  are  not  pertinent  to  our  analysis,  but  may  be  found  in  the 
listed  references;  the  pertinent  aspects  are  provided  below. 

1.  Bathymetry 

Figure  7  shows  the  bathymetry  of  the  acoustic  channel,  as  well  as  the 
locations  of  the  source  (far  left  side  of  the  plot)  and  the  receiver  (located  at 
roughly  34  km  range)  (to  be  discussed  later).  The  geometric  axis  from  the 
source  to  the  receiver  was  almost  directly  downslope. 

2.  Sound  Speed  Profile 

Figure  7  and  Figure  8  illustrate  the  sound  speed  field  in  the  channel. 
The  curves  in  Figme  8  were  obtained  by  interpolation  of  sensor  casts  made  at 
roughly  10  km  intervals;  the  dotted  and  solid  cmwes  correspond  to  the  source 
and  receiver  locations,  respectively.  The  sound  speed  field  obtained  from  the 
interpolation  was  used  as  input  to  the  BBCM  model  for  all  data  analysis.  The 
higher  resolution  sound  speed  field  of  Figm*e  7  was  obtained  by  tomographic 
inversion  [Ref.  19];  it  is  provided  for  illustration  only  and  was  not  used  in  our 
analysis. 
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Figure  7:  Sound  speed  field  and  bathymetry 


The  plot  uses  MATLAB’s  “pseudocolor”  feature:  the  gray  level  at  each 
point  in  the  plot  maps  to  a  sound  speed  (in  meters  per  second)  per  the  gray 
level  bar  at  the  left  side  of  the  plot.  The  plot  clearly  shows  the  front  which 
was  the  primary  subject  of  interest  in  the  Barents  Sea  Polar  Front 
Experiment;  this  front  was  nearly  perpendicular  to  the  axis  between  the 
source  and  receiver.  This  fact,  combined  with  the  fact  that  sound  propagation 
was  almost  entirely  downslope,  allows  us  to  make  the  assumption  that  no 
horizontal  refraction  or  azimuthal  scattering  occurred  [Ref.  19].  The  front 
was  observed  to  move  upslope  and  downslope  with  a  dominant  periodicity  of 
about  two  hours  and  a  peak-to-peak  amplitude  on  the  order  of  4  km.  The 
sound  speed  and  density  in  the  bottom  were  obtained  from  standard  Navy 
databases  and  were  found  to  be  3200  m/s  and  2600  kg/m^,  respectively;  these 
values  were  verified  by  SUS  measurements  [Ref.  53]. 
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Figure  8:  Sound  Speed  Profiles 


B.  SIGNAL  AND  NOISE 

The  transmitted  signal  consisted  of  M-sequences  with  a  center 
frequency  of  224  Hz.  Two  sets  of  M-sequences,  separated  by  about  nine  hours, 
were  transmitted  during  the  experiment.  Each  set  consisted  of  30  M- 
sequences,  each  of  about  five-second  duration,  giving  a  total  of  about  2.5 
minutes  of  transmission  per  set.  The  unique  properties  of  the  M-sequence 
were  needed  to  achieve  the  goals  of  the  Barents  Sea  Experiment  (i.e., 
accurate  travel  time  determination),  but  are  not  relevant  to  this  dissertation 
and  will  therefore  not  be  discussed  here.  For  the  present  analysis,  the 
received  signal  was  filtered  using  a  Minimum- Variance  filter  to  remove  the 
M-sequence  properties;  this  type  of  filter  was  selected  to  ensure  that  any 
strong  interfering  signals  due  to  engine  noise  from  the  test  ship  were  nulled 
out  (see,  e.g.,  [Ref.  54]).  We  found,  not  surprisingly,  that  use  of  this  filter 
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resulted  in  emitter  location  estimates  superior  to  those  obtained  using  the 
usual  Fast  Fourier  Transform  (FFT)  approach.  Figure  9  shows  the  power 
spectral  density  (relative  to  the  peak)  of  the  unfiltered  received  signal  at  a 
particular  hydrophone  and  a  particular  time.  The  received  signal  (at  224  Hz), 
when  averaged  over  all  hydrophones,  has  a  SNR  of  about  10  dB  (ratio  taken 
over  the  entire  signal  bandwidth);  the  exact  value  of  the  SNR  is  not 
particularly  important  for  our  analysis. 


Signal  Power  Spectrum 


Figure  9:  Power  Spectrum 


C.  RECEIVER 

The  receive  array  consisted  of  a  vertical  string  of  16  identical, 
omnidirectional  hydrophones  with  10  m  spacing.  The  uppermost  hydrophone 
was  at  a  depth  of  123.8  m.  Prior  to  the  collection  of  data  used  in  this 
dissertation,  hydrophones  1,  2,  4,  6-8,  and  16  experienced  partial  failure 
(separation  of  the  two  piezoelectric  cylinders  comprising  each  hydrophone) 
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which  reduced  their  sensitivity  by  6  dB,  This  sensitivity  reduction  was 
incorporated  into  the  data  analysis. 
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IV.  ANALYSIS  AND  RESULTS 


This  chapter  provides  further  details  of  the  analytical  technique, 
including  preprocessing  methods.  It  also  provides  and  interprets  our  main 
results  (which  take  the  form  of  ambiguity  surface  plots)  from  application  of 
the  MFP/MMP  estimators  to  the  Barents  Sea  data  for  various  choices  of 
parameters  (noise,  data  length,  etc.). 

A  ANALYSIS  TECHNIQUE 

1.  Signal  Processing 

The  computational  procedure  may  be  outlined  as  follows: 

•  Construct  the  density  and  sound  speed  fields  (functions  of  range 
and  depth)  from  in  situ  measurements  using  bilinear  interpolation); 

•  Using  sound  speed,  density,  receiver  horizontal  location,  and 
bathymetry  as  inputs,  and  regarding  each  receive  hydrophone  as  a 
unit  emitter,  generate  (using  a  suitable  propagation  model)  and 
store  the  parameters  appearing  in  (31)  (a  separate  set  of 
parameters  for  each  hydrophone  depth); 

•  Using  a  numerical  implementation  of  (31)  and  invoking  the 
principle  of  acoustic  reciprocity,  calculate  the  array  manifold  on  a 
suitable  grid  of  (royZo)  values; 

•  Take  the  raw  hydrophone  data  (an  M-sequence)  and  add  the 
desired  amount  of  noise  (measured  in  the  same  environment  during 
a  period  when  no  transmissions  occurred); 
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•  Pass  the  resulting  signal  through  a  50th  order  minimum-variance 
filter  to  eliminate  out-of-band  noise  and  to  remove  the  M-sequence 
properties  from  the  signal  (i.e.  make  it  into  an  ordinary  narrow- 
band  process); 

•  Calculate  estimated  spatial  and  modal  covariance  and  cumulant 
matrices  from  the  resulting  signal; 

•  Match  these  covariance  matrices  against  the  replica  fields 
calculated  above  (for  MMP,  use  only  the  components  corresponding 
to  the  desired  modes). 

All  computations  were  done  using  Matlab  4.2c  on  a  Hewlett-Packard 
735  Workstation.  The  propagation  modeling  used  the  Broadband  Coupled 
Mode  (BBCM)  algorithm  [Ref.  19].  In  every  case,  the  ambiguity  surfaces  were 
calculated  for  a  grid  with  the  following  specifications,  selected  to  ensure 
peaks  would  not  be  missed  while  keeping  the  computational  load  reasonable: 


Parameter 

Minimum 

Maximum 

Increment 

(Emitter)  Range 

15  km 

40  km 

40  m 

Depth 

2m 

146  m 

2m 

Table  1:  Computation  grid 

The  emitter  range  used  in  the  plots  is  measured  with  respect  to  a  reference 
different  from  those  shown  in  Figure  6  and  Figure  7:  it  is  measured  with 
respect  to  a  point  1675  m  downslope  from  the  receiver.  The  15^0  km  range 
window  was  selected  to  allow  a  fair  assessment  of  our  estimators,  while 
keeping  the  amount  of  computation  manageable.  Because  the  water  gets 
deeper  as  one  gets  closer  to  the  receiver,  more  modes  are  required  to 
construct  the  pressure  field  (the  number  of  modes  supported  is  roughly 
proportional  to  the  water  depth  for  a  fixed  frequency).  The  run  time  required 
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by  the  BBCM  model  appears  to  be  proportional  to  the  number  of  modes 
cubed. 

2.  Plotting 

The  ambiguity  surface  plots  were  generated  using  Matlab  (version 
4.2a)  on  a  Macintosh  LCIII  and  printed  on  a  600  dot  per  inch  (dpi)  HP 
LaserJet  4,  using  the  Matlab  “pseudocolor”  plot  function  with  32  gray  levels. 
Because  a  600  dpi  printer  is  not  able  to  generate  a  dot  screen  with  adequate 
resolution  to  display  data  on  a  grid  as  large  as  that  found  in  Table  1  without 
requiring  an  excessive  amount  of  space  on  the  page,  the  ambiguity  surface 
data  was  smoothed  before  plotting:  at  each  depth  grid  point,  the  value  plotted 
is  the  average  for  three  adjacent  range  grid  points.  Subjectively  speaking, 
little  information  appeared  to  be  lost  as  a  result  of  this  smoothing. 

The  quantitative  assessment  of  the  effect  of  various  parameters  on 
localization  performance  presents  an  interesting  problem.  On  all  plots,  the 
bright  areas  correspond  to  maxima  of  the  functions  (27)  through  (30).  No  gray 
scale  is  provided,  since,  for  the  MUSIC  algorithm,  the  peak  height  does  not 
necessarily  correspond  to  the  power  of  the  received  signal.  In  fact,  the  values 
corresponding  to  “white”  and  “black”  differ  somewhat  from  plot  to  plot.  A 
more  suitable  approach  for  quantitatively  comparing  the  different  plots  may 
be  based  on  three  major  performance  measures  used  in  DOA  estimation 
literature  (see,  e.g.,  [Ref  23]): 

•  Bias  in  the  emitter  location  estimates; 

•  Ability  to  resolve  multiple  closely  spaced  emitters;  and 

•  Existence  of  peaks  at  locations  where  no  emitters  exist. 
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In  the  present  case,  because  the  precise  location  of  the  emitter  is  not  known, 
the  first  measiu-e  is  not  applicable  in  the  form  stated.  Also,  only  one  emitter 
is  present,  so  the  second  measure  is  not  particularly  useful  either.  The  third 
measure  is  applicable,  although  not  easy  to  quantify.  We  have  elected  instead 
to  use  the  following: 

•  Ratio  of  the  height  of  the  correct  peak  to  that  of  the  highest  false 
peak  (Ml); 

•  Size  (area)  of  the  correct  peak  relative  to  the  area  of  the  entire  grid 
(M2); 

•  Ratio  of  the  average  height  of  all  false  peaks  to  the  height  of  the 
correct  peak  (M3). 

Obviously,  when  M1<1,  an  incorrect  estimate  for  the  emitter  location  will  be 
obtained.  The  values  of  these  three  measures  are  given  in  the  caption  for 
each  figure  that  follows. 

B.  OVERVIEW  OF  RESULTS 

In  each  section  to  follow,  we  consider  the  effect  of  the  variation  of  a 
single  parameter  {e.g.,  data  length,  algorithm  t3rpe,  SNR,  etc.)  while  holding 
all  other  parameters  fixed  at  the  following  nominal  values  (for  which 
performance  is  good): 

•  Matched-mode  method  (using  projection  method  of  mode  filtering 
and  modes  1-4); 

•  Data  length  of  1024  points; 
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•  Data  taken  from  the  beginning  of  the  second  set  of  M-sequence 
transmissions; 

•  SNR=10  dB  (i.e.,  no  noise  beyond  that  actually  observed  with  the 
signal); 

•  MUSIC  method  of  array  processing; 

•  Second-order  statistics; 

•  Coupled-mode  propagation  model;  and 

•  Hydrophone  data  corrected  for  reduced  sensitivity  of  damaged 
phones. 

Each  plot  to  follow  has  a  title  at  the  top  containing  the  most 
significant  parameters  (MFP/MMP,  SNR,  data  length,  and  array  processing 
algorithm).  Additional  pertinent  information  appears  either  in  the  caption  or 
in  the  text  referring  to  the  plot.  In  every  case,  the  correct  emitter  location  is 
at  approximately  36  km  range  and  122  m  depth. 

C.  MATCHED-FIELD  VERSUS  MATCHED-MODE  PROCESSING 

As  mentioned  earlier,  propagation  models  are  generally  more  sensitive 
to  errors  in  knowledge  of  the  environment  when  predicting  high-order  modes 
than  low-order  modes.  We  therefore  expect  MFP  to  exhibit  a  higher  incidence 
of  false  peaks  than  MMP,  since  MFP  uses  all  available  modes  (26).  Figure  10 
(MFP)  and  Figure  11  (MMP)  illustrate  this  effect.  For  purposes  of 
illustration,  we  have  added  a  small  white  circle  at  the  correct  emitter  location 
in  Figure  10,  although  we  have  not  done  so  with  subsequent  plots,  because 
the  emitter  location  is  fixed.  Although  MFP  gives  a  peak  at  approximately 
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the  correct  emitter  location  (Ml=0.94),  there  are  2  other  peaks  which  are 
higher,  as  well  as  numerous  smaller  peaks  (M3=0.043).  MMP  gives  the 
highest  peak  at  the  correct  location  (Ml=1.7)  and  has  fewer  and  smaller  false 
peaks  (M3=0.016). 

MF,  MUSIC,  SNR=+10  dB,  1024  pts 
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Figure  10:  Matched-Field  Processing  (Ml=0.94;  M2=0.0077;  M3=0.043) 
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MM,  MUSIC,  Modes  1-4,  SNR=+10  dB,  1024  pts 
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Figure  11:  Matched-Mode  Processing  (Ml=1.7;  M2s0.012;  M3=0.016) 

D.  BEHAVIOR  OF  DIFFERENT  ARRAY  PROCESSING  METHODS 

Both  Figure  12  (MVM)  and  Figure  13  (Bartlett)  show  poor  resolution 
compared  with  MUSIC  (Figure  11).  Both  of  these  methods  produce  the 
largest  peak  at  the  correct  location  (Ml=1.26  and  1.03,  respectively),  but 
there  are  large  and  numerous  false  peaks  (M3=0.11  and  0.14,  respectively). 
In  particular,  with  the  Bartlett  method,  the  correct  peak  is  nearly  impossible 
to  identify  visually.  The  behavior  of  these  three  methods  when  used  on  this 
data  set  is  thus  consistent  with  that  observed  in  DOA  estimation  and  with 
modeled  MFP  data  (see,  e.g.,  [Refs.  3,  23]). 


f: 
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Depth  (m)  Depth  (m) 


MM,  MVM,  Modes  1-4,  SNR=+10  dB,  1024  pts 


Range  (m)  ^ 

12;  Minimum-Variance  Method  (Ml=1.26;  M2=0.033;  M3=0.11) 


MM,  Bartlett,  Modes  1-4,  SNR=+10  dB,  1024  pts 
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Figure  13;  Bartlett  Method  (Ml=1.03;  M2=0.017;  M3=0.14) 
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E.  EFFECT  OF  DATA  LENGTH  ON  PERFORMANCE 


As  the  length  of  data  increases  (assuming  temporally  stationary  data 
and  fixed  SNR),  the  accuracy  of  the  estimate  of  spatial  covariance  should 
improve.  On  the  other  hand,  if  significant  non-stationarity  is  present,  we 
expect  that  increases  in  data  length  past  a  certain  point  (depending  on  SNR) 
may  actually  degrade  the  accuracy  of  this  estimate.  This  is  the  case  with  our 
data,  as  may  be  seen  in  Figure  14  through  Figure  16,  where  the  ambiguity 
surface  for  a  data  length  of  512  shows  some  improvement  (Ml=2.2, 
M3=0.015)  with  respect  to  the  surface  for  a  data  length  of  1024  (Figure  11), 
and  the  size  of  the  false  peaks  increases  noticeably  for  data  lengths  of  2048 
(Ml=1.6,  M3=0.019)  and  4096  (Ml=1.4,  M3=.020).  This  behavior  is 
presumably  due  to  surface  wave  effects,  which  are  the  only  likely  source  of 
temporal  variability  over  the  roughly  one-second  time  scale  involved  here. 
Obviously,  data  length  can  only  be  reduced  so  far  before  the  benefits  gained 
by  avoiding  non-stationarity  are  outweighed  by  estimation  errors  arising 
fi*om  low  information-to-noise  ratio.  In  fact,  performance  degradation  became 
evident  at  SNR=+10  dB  when  the  data  length  was  reduced  to  256  points  (not 
shown). 
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MM,  MUSIC,  Modes  1-4,  SNR=+10  dB,  512  pts 
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Figure  14:  Data  Length  512  (Ml=2.2;  M2=0.012;  M3=0.015) 
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Figure  15:  Data  Length  2048  (Ml=1.6;  M2=0.012;  M3=0.019) 
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MM.  MUSIC,  Modes  1-4,  SNR=+10  dB,  4096  pts 
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Figure  16:  Data  Length  4096  (Ml=1.4;  M2=:0.012;  M3a0.020) 

F.  EFFECT  OF  NOISE  ON  PERFORMANCE 

Figure  17  through  Figure  19  illustrate  the  effect  of  additive  noise. 
These  plots  exhibit  high  false  peaks  compared  with  Figure  11  (no  added 
noise).  Below  a  SNR  of  —10  dB,  the  MUSIC  method  using  2nd  order  statistics 
no  longer  gives  the  largest  peak  at  the  actual  emitter  location  (Ml=0.72  for 
Figure  18).  Interestingly,  the  MUSIC  method  with  4th-order  statistics 
actually  gives  poorer  results  than  with  2nd-order  statistics  at  all  SNRs 
(Figure  19  shows  the  0  dB  result).  One  reason  for  this  somewhat  surprising 
result  appears  to  be  that  the  signal  turns  out  to  be  roughly  as  close  to 
Gaussian  as  the  additive  noise,  as  measured  by  the  difference  between  its 
4th-order  moment  and  the  4th-order  moment  of  a  Gaussian  process  with  the 
same  lower-order  moments  (11).  Specifically,  let  us  define  a  measure  G  for 
quantifying  Gaussianity  for  the  received  signal  as  follows: 
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J!£iL 


(40) 


where  ||- 1|^  indicates  the  Frobenius  norm*,  C4  is  defined  by  (14),  and  is  the 
matrix  of  4th  order  moments;  i.e., 

M.^{i,j)  =  E[piP*PjP*]. 


As  discussed  earlier,  if  the  processes  are  Gaussian,  €4=0  and  therefore 
G=0.  We  find  that,  in  this  experiment,  for  noise  alone,  G=0.13  and  for  signal 
alone  G=0.12  (the  use  of  different  norms  in  (40)  does  not  significantly  affect 
these  values).  Recall  that  our  motivation  for  using  higher-order  statistics  in 
the  first  place  was  based  on  an  expectation  that  the  noise  would  be  Gaussian 
and  the  signal  would  not.  It  appears  that  this  Gaussian  property  is 
characteristic  of  the  filtered  M-sequence  signal;  for  continuous  wave  (CW) 
data  gathered  later  in  the  experiment  (at  which  time,  unfortunately, 
environmental  measurements  are  not  available),  G  is  significantly  higher. 


*  The  Frobenius  norm  of  a  matrix  A  is  defined  as  ||A||^ 
components  of  A. 


si4r  ,  where  Ay  are  the 
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MM,  MUSIC,  Modes  1-4,  SNR=:+0  dB,  1024  pts 
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Figure  19:  SNR  =  0  dB,  4th  order  statistics 
(Ml=0.34;  M2=0.008;  M3=0.015) 

A  more  significant  flaw  in  the  cumulant-based  MUSIC  estimator  is 
illustrated  in  Figure  20  through  Figure  23,  which  show  the  performance  of 
cumulant-based  MUSIC  MFP  versus  that  of  conventional  MUSIC  MFP  when 
applied  to  a  synthetic  data  set  (no  mismatch  between  the  actual  and  predicted 
sound  fields)  for  SNRs  of  +10  dB  and  0  dB.  At  +10  dB,  the  cumulant-based 
method  (Ml=12.3)  greatly  outperforms  the  conventional  method  (Ml=2.9). 
However,  at  0  dB,  the  conventional  method  is  superior  (Ml=2.7  vs.  Ml=1.6), 
despite  the  fact  that  the  signal  subspace  is  estimated  much  more  accurately 
for  the  cumulant  method  than  for  the  conventional  method  (as  quantified  by 
the  angle  between  the  signal  subspaces  at  +10  dB  and  0  dB).  This  behavior  is 
due  to  the  fact  that  both  the  cumulant  matrix  C4  and  the  array  manifold 
vectors  (aoa*)  are  defined  to  be  real:  since  both  MFP  and  MMP  rely  heavily  on 
phase  information  for  accurate  localization,  use  of  amplitude  information 
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alone  causes  the  estimators  to  be  highly  sensitive  to  noise-induced  errors  in 
the  estimate  of  C4. 


MF,  MUSIC,  SNR=+10  dB,  1024  pts 


Figure  20:  Conventional  MUSIC,  synthetic  data 
(Ml=2.9,  M2=0.0091,  M3=0.014) 
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Depth  (m)  Depth  (m) 


MF,  MUSIC,  SNR=+10  dB,  1024  pts 
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Figure  21:  Cumulant  MUSIC,  synthetic  data 
(Ml=12.3,  M2=0.0021,  M3=0.0022) 


MF,  MUSIC,  SNR=+0  dB,  1024  pts 


Range  (m)  ^  ^0" 

Figure  22:  Conventional  MUSIC,  synthetic  data 
(Ml=2.7,  M2=0.0093,  M3=0.015) 
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MF,  MUSIC,  SNR=+0  dB,  1024  pts 


Figure  23:  Cvunulant  MUSIC,  synthetic  data 
(Ml=1.6,  M2=0.0022,  M3=0.015) 

G.  BEHAVIOR  OF  DIFFERENT  MODE  INVERSION  METHODS 

Figure  24  shows  the  result  when  the  pseudoinverse  method  is  used 
instead  of  the  projection  method.  As  mentioned  above,  the  two  methods  are 
very  similar  mathematically  and  give  about  the  same  performance  (compare 
with  Figure  11).  The  false  peaks  are  slightly  larger  and  more  numerous  with 
the  pseudoinverse  method  (M3=0.020  vs.  0.016  for  the  projection  method). 
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MM,  MUSIC,  Modes  1-4,  SNR=+10  dB,  1024  pts 
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Figure  24:  Pseudo-inverse  mode  filter  (Ml=1.7;  M2=0.013;  M3=0.020) 

H.  EFFECT  OF  ARRAY  SHADING 

Figure  25  shows  the  effect  of  not  including  the  required  sensitivity 
correction  for  the  failed  hydrophones  (i.e.,  no  array  shading).  The  higher  false 
peaks  are  appeirent  (Ml=l.l  and  M3=0.025,  as  compared  with  Ml=l,7  and 
M3=0.016  when  shading  is  used). 
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MM,  MUSIC,  Modes  1-4,  SNR=+10  dB,  1024  pts 


Range  (m) 


Figure  25:  No  phone  sensitivity  correction 
(Ml=l.l;  M2=0.012;  M3=0.025) 

I.  TEMPORAL  VARIABILITY  OF  RESULTS 

Figure  26  shows  the  results  from  a  data  segment  obtained  during  the 
first  set  of  M-sequence  transmissions.  Although  there  is  a  peak  at  the  correct 
emitter  location,  it  is  not  the  largest  peak  (Ml=0.46).  The  change  in 
localization  performance  with  respect  to  that  obtained  with  a  data  segment 
from  the  second  set  of  transmissions  is  probably  due  to  temporal  fluctuations 
in  the  sound  speed  field,  since  none  of  the  other  physical  parameters  changed 
significantly  between  the  two  data  sets.  This  behavior  is  not  surprising,  since 
the  period  of  frontal  motion  (two  hours)  is  much  less  than  the  time  between 
the  sets  of  M-sequence  transmissions.  Although  the  plots  are  not  shown  here, 
we  found  that  the  localization  performance  remained  qualitatively  the  same 
for  all  data  segments  within  a  given  set  of  transmissions. 
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Figure  26:  Data  from  1st  transmission  (Ml=0.46;  M2=0.015;  M3=0.026) 

J.  EFFECT  OF  MODE  SELECTION  ON  MMP  RESULTS 

As  discussed  above,  it  is  important  to  choose  a  suitable  mode  subset  to 
ensure  satisfactory  results.  Obviously,  with  60  modes  available,  there  is  a 
very  large  number  of  potential  combinations,  only  a  few  of  which  will  be 
presented  here.  Figure  27  shows  the  result  when  only  the  first  3  modes  are 
used;  a  small  peak  is  visible  at  the  correct  location  (Ml=0.46).  Figure  28 
shows  the  result  when  modes  1-5  are  used;  the  plot  shows  a  slight 
improvement  in  performance  compared  to  the  result  with  modes  1-4  (Ml=1.8 
and  M3=0.015  vs.  1.7  and  0.016,  respectively).  Using  more  than  the  first  five 
modes  tends  to  increase  the  number  and  size  of  the  false  peaks.  Figure  29 
(modes  1-6  used)  shows  this  effect  (Ml=1.3,  M3=0.018). 
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.5  2  2.5  3  3.5 

Range  (m)  x  10“ 


68 


MM,  MUSIC,  Modes  1-6,  SNR=+10  dB,  1024  pts 
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Figure  29:  Modes  1-6  (Ml=1.3;  M2=0.011;  M3=0.018) 


K  EFFECT  OF  MODEL  SELECTION  ON  PERFORMANCE 

Figure  30  shows  the  result  when  the  mode  coupling  accounted  for  by 
the  BBCM  method  is  not  included  in  calculation  of  the  replica  fields,  that  is, 
when  the  range  dependence  of  the  replica  fields  arises  only  from  the  range 
dependence  of  the  mode  functions  and  wavenumbers  (the  adiabatic 
approximation  of  (22)).  There  is  no  indication  of  a  peak  at  or  near  the  actual 
emitter  location,  so  the  performance  measures  Ml,  M2,  M3  are  meaningless 
and  are  not  provided. 
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Figure  30:  Adiabatic  model 
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V.  CONCLUSIONS 


The  results  presented  in  the  body  of  this  dissertation  clearly 
demonstrate  that  MUSIC-based  MMP  techniques  may  be  effectively 
employed  even  in  a  very  challenging  acoustic  environment  such  as  that  which 
existed  during  the  Barents  Sea  Polar  Front  Experiment.  It  is  appropriate  to 
review  some  of  the  unique  features  of  this  environment  (as  compared  to 
idealized  numerical  simulations  or  simple  experiments): 

•  Strong  range  dependence  of  bathymetry  and  the  sound  speed  field 
(which  included  a  strong,  rapidly  moving  front); 

•  A  degraded  receive  array  spanning  about  half  the  water  column 
and  having  many  fewer  elements  than  the  number  of  modes 
supported  by  the  channel; 

•  Relatively  coarse  sampling  of  the  sound  speed  field  (only  about  once 
per  10  km  interval,  notwithstanding  the  presence  of  a  front); 

Despite  these  challenges,  we  achieved  considerable  success  with  our 
localization  approach.  Some  of  the  specific  findings  and  original  contributions 
of  this  research  are: 

•  Contrary  to  much  conventional  wisdom,  the  subspace-based  MUSIC 
method  produced  good  results  despite  the  inaccuracies  inherent  in 
experimental  data.  In  fact,  the  MUSIC  algorithm’s  high  resolution 
was  vital  to  accurate  localization,  because  the  receive  array  was 
relatively  short  and  few  robust  modes  were  available. 
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•  The  Broad  Band  Coupled  Mode  (BBCM)  model  generates  replica 
fields  with  sufficient  accuracy  to  allow  localization  via  MMP  in  this 
strongly  range-dependent  environment.  The  adiabatic  approxima¬ 
tion  was  found  to  be  grossly  inadequate  for  this  environment. 

•  The  cumulant-based  MUSIC  estimator  used  in  this  dissertation 
was  too  sensitive  to  noise-induced  and  model-induced  estimation 
errors  to  be  useful  with  real  data.  This  behavior  was  due  to  the  fact 
that  the  cumulant  matrix  and  the  replica  fields  were  defined  as  real 
quantities;  thus,  the  phase  information,  which  is  vital  to  robust  and 
accurate  localization,  was  not  available. 

•  In  an  environment  with  strong  temporal  variability,  localization 
performance  can  vary  drastically  over  relatively  short  time  scales. 

In  summary,  our  approach,  which  combined  the  high-resolution  MUSIC 
algorithm  with  MMP,  allowed  accurate  localization  estimates,  even  though 
only  a  few  robust  modes  could  be  obtained  via  mode  filtering. 

The  approach  used  in  our  analysis  may  be  modified  in  three  obvious 
respects,  each  of  which  offers  potentially  significant  improvement  in 
localization  performance  and  is  worthy  of  further  study. 

The  first  modification  relates  to  the  assumptions  concerning 
observation  noise.  The  basic  MUSIC  algorithm  used  in  our  research  assumes 
that  the  noise  covariance  is  some  multiple  of  the  identity  matrix  (i.e., 
spatially  isotropic).  As  mentioned  earlier,  the  MUSIC  algorithm  may  be 
extended  to  the  case  where  the  noise  is  not  isotropic  via  use  of  a  generalized 
eigendecomposition.  Use  of  this  extended  MUSIC  algorithm  has  the  potential 
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to  further  lower  the  SNR  threshold  (the  SNR  above  which  the  estimator 
performs  satisfactorily)  observed  in  our  results. 

A  second  way  to  improve  localization  performance  involves 
modification  of  our  cumulant-based  MUSIC  estimator.  We  noted  earlier  that 
an  extended  version  of  the  cumulant  matrix  has  been  defined  [Ref.  20]. 
Although  an  estimator  based  on  this  matrix  would  be  more  computationally 
intensive  than  the  estimator  defined  here,  it  is  expected  to  be  less  sensitive  to 
estimation  errors. 

Refinement  in  the  process  of  selecting  suitable  modes  for  MMP 
localization  is  a  third  means  of  improving  the  estimator.  Although  the  simple 
approach  described  here  was  effective  in  generating  an  appropriate  mode  set 
for  this  environment,  it  may  not  produce  results  of  the  same  quality  in  other 
environments.  An  approach  rel5dng  more  heavily  on  the  propagation  physics 
of  the  channel  could  greatly  reduce  the  amount  of  “trial  and  error”  involved  in 
the  process. 
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